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This paper focuses on unconditionally optimal error analysis of an uncoupled and 
. linearized Crank-Nicolson Galerkin finite element method for the time-dependent non- 

linear thermistor equations in <i-dimensional space, d = 2, 3. We split the error function 
into two parts, one from the spatial discretization and one from the temporal dis- 
i-pH " cretization, by introducing a corresponding time-discrete (elliptic) system. We present 

a rigorous analysis for the regularity of the solution of the time-discrete system and er- 
ror estimates of the time discretization. With these estimates and the proved regularity, 
optimal error estimates of the fully discrete Crank-Nicolson Galerkin method are ob- 
tained unconditionally. Numerical results confirm our analysis and show the efficiency 
. of the method. 
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1 Introduction 

We consider the time-dependent nonlinear thermistor system 



^-A U = a(u)\Vd ) \ 2 , (1.1) 



- V • (<t(u)V0) = 0, (1.2) 

for x G 17 and t € [0,T], where Q is a bounded domain in M d , d = 2,3. The initial and 
boundary conditions are given by 

u(x,t)=0, 4>(x,t) = g(x,t) for x £ <9f2, t e [0,T], 
u(x,0) = uo(x) for x G Q. 

The nonlinear system above describes the model of electric heating of a conducting body, 
where u is the temperature, (f> is the electric potential, and a is the temperature-dependent 
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electric conductivity. Following the previous works 138] . we assume that a £ W 1,00 (M) 
and 

ai < a(s) < a 2 , (1.4) 

for some positive constants <7i and 02- 

Theoretical analysis for the time-dependent thermistor equations was done by several 
authors [3J [TOj [361 EZj- Among these works, Yuan and Liu [37] proved the existence 
and uniqueness of a C a solution in three-dimensional space. Based on their result, further 
regularity can be derived with suitable assumptions on the initial and boundary conditions. 
Numerical methods and analysis for the thermistor system can be found in [2J [H [TU [351 ESI 
139] , For the system in two-dimensional space, the optimal L 2 error estimate of a mixed finite 
element method with a linearized semi-implicit Euler scheme was obtained in |38] under a 
weak time-step condition. Error analysis for the three-dimensional model was given in [14] , 
in which a linearized semi-implicit Euler scheme with a linear Galerkin FEM was used. An 
optimal L 2 -error estimate was obtained under the condition r = 0{h 1 / 2 ). A more general 
time discretization with higher-order finite element approximations was studied in [2]. An 
optimal L 2 -norm error estimate was given under the conditions r = 0(h 3 ^ 2p ) and r > 2, 
where p is the order of the time discretization and r is the degree of piecewise polynomials 
of the finite element space. 

Clearly, there are several different time discretizations for nonlinear parabolic systems, 
explicit, semi-explicit (or semi-implicit) and implicit. The most popular and widely-used 
approach is linearized (semi)-implicit scheme. At each time step, the scheme only requires 
the solution of a linear system. However, time-step condition is always a key issue for such a 
scheme. To study the error estimate of linearized (semi)-implicit schemes, the boundedness 
of the numerical solution (or error function) in L°° norm or a stronger norm is often required. 
If a priori estimate for numerical solution in such a norm cannot be provided, one may 
employ the mathematical induction with an inverse inequality to bound the numerical 
solution, such as, by 

\\U% - R h u n \\ L „ < Ch- d l 2 \\Ul - R h u n \ L 2 < Ch- d ' 2 {T m + h r+1 ), (1.5) 

where UJf is the finite element solution, u n is the exact solution and Rh is certain projection 
operator. The above approach, however, requires a time-step condition r = 0(h d ^ 2m ^). 
This approach has been widely used in the error analysis of many different nonlinear 
parabolic PDEs, e.g., see [161 [13 [23] for Navier-Stokes equations, [21 Q31 [38] for nonlin- 
ear thermistor problems, [151 [281 ETJ for porous media flows, [HI [32] for viscoelastic fluid 
flow, |24] for KdV equations, [91 [25] for the Ginzburg-Landau equations, [61 [30] for nonlin- 
ear Schrodinger equations and [12} [34] for some other equations. In all these works, error 
estimates were established under certain time step restrictions. The time-step restrictions 
arising from theoretical analysis may result in the use of a very small time step and ex- 
tremely time-consuming in practical computations. However, we believe that such time-step 
conditions may not be necessary for most cases. A new approach was introduced in our 
recent works [20 \ [2T] , also see [22] , in which the error estimates of a linearized backward 
Euler Galerkin methods for a porous media flow and the thermistor system were obtained, 
respectively, under the condition of h and r being smaller than a positive constant. In 
this paper, we propose an uncoupled and linearized Crank-Nicolson Galerkin finite ele- 
ment method for the nonlinear thermistor system and present optimal error estimates in 
both L 2 and H 1 norms without any stepsize restrictions. In this method, the standard 



2 



Crank-Nicolson scheme is applied for the linear term in the temperature equation and an 
extrapolation approximation is used for the nonlinear electric conductivity. At each time 
step, one only needs to solve two uncoupled linear systems. The main idea of our aprooach 
is to split the error function into two parts, the spatially discrete error and the temporally 
discrete error, by introducing a corresponding time-discrete (elliptic) system. The former 
arises from the Galerkin FEM discretization for the time-discrete equations and depends 
only upon the spatial mesh size h (independent of the time-step size r). If a suitable regu- 
larity of the solution to the time-discrete equations has been proved, the numerical solution 
can be bounded by 

\\U% - R h U n \\ L oo < Ch- d l 2 \\Ul - R h U n \\ L 2 < Ch~ d l 2 h r+1 (1.6) 

without any time-step condition, where U n is the solution of the time-discrete equations. 
More important is that our approach is applicable for more general nonlinear parabolic 
PDEs and many other time discretizations to obtain unconditional convergence and optimal 
error estimates. 

The rest of the paper is organized as follows. In Section 2, we present the uncoupled and 
linearized Crank-Nicolson scheme with a linear Galerkin finite element approximation in 
the spatial direction and state our main results. After introducing the corresponding time- 
discrete system, we provide in Section 3 a priori estimates and optimal error estimates for 
the time-discrete solution, which imply the suitable regularity of the time-discrete solution. 
With the regularity obtained, in Section 4 we present optimal error estimates of the fully 
discrete Galerkin finite element solution in both the L 2 norm and the H l norm without any 
time-step conditions. Numerical results are presented in Section 5 to confirm our theoretical 
analysis. 

2 The main result 

Let Q be a bounded, smooth and convex domain in M. d (d = 2,3). Let 7T/j be a regular 
division of ft into triangles Tj, j = l,--- , M in K 2 or tetrahedra in ]R 3 , and denote by 
h = maxi<j<jvf{diamTj} the mesh size. For a triangle (or tetrahedra) Tj at the boundary, 
we define Tj to be a triangle with one curved side (or a tetrahedra with one curved face 
in R 3 ) with the same vertices as Tj, and set Dj = Tj\Tj. For an interior triangle, we set 
Tj = Tj and Dj = 0. For a given triangular (or tetrahedral) division of fi, we define the 
finite element spaces [29J: 

Vh = {vh £ C(Q) : Vh\Tj is linear at each element and Vh = on Dj}, 
Sh = {vh £ C(O) : Vh\f. is linear at each element}. 

It follows that Vh is a subspace of Hq(Q) and is a subspace of i? 1 (0). For any function 
v G Sh, we define A/jf to be a function satisfying K^v = on Dj and AhV = v on Tj. 
We further define 11^ : C${Q) — > Sh to be the Lagrangian interpolation operator and set 
lih = A/jII^. Clearly, IT^ is a projection operator from Co(O) onto Vh- 

Let = to < *i < • • • < tN = T be a uniform partition of the time interval [0,T] with 
t n = tit and let 

u n = u(x,t n ), <p n = (j)(x, t n ) for n = 0, 1, • • • ,N. (2.1) 
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For any sequence of functions {/ n }„ =0 ' we define 

D T f n+l = J - *— J n+1 ' 2 = (3/ m - r -1 )/2, (2.2) 

r 

P +1/2 = l(f n + f n+1 ), (2.3) 

for n = 1,2, - •• ,JV-1. 

For the simplicity of notations, we denote by C a generic positive constant and by e a 
generic small positive constant, which depend solely upon the physical parameters of the 
problem and independent of t, h and n. We assume that g(-,t) G H 1 ^) is given for each 
fixed t > 0. 

We propose an uncoupled and linearized Crank-Nicolson Galerkin finite element method 
to solve the system ([13]) - ([TH]) . which seeks U™ +1 G V h and $^ +1/2 G £f n ~ 1/2 + 14, n = 
0, 1, • • • , iV - 1, such that 



a{U n h +1/2 )V<5>l +1/ \ ) = 0, V^y, (2.4) 



D T ltf+\«) + (W^ 2 , W) = ( CT (C/; +1/2 )|V^ +1 / 2 | 2 ,,;), V « G V h (2.5) 



where a standard extrapolation [13] is used to approximate the nonlinear electric conduc- 
tivity for n > 0. 

At the initial time steps, we choose U® = Il^ito and let be the Galerkin solution to 
the potential equation 

v(uo)V$° h ,V(p) =0, V y G 14- (2.6) 



^1/2 

and can be calculated either by a semi-implict Euler scheme 
Pr 1 / 2 



(^-— — , U ) + (vt/; /2 ) V U ) = ( < 7( U0 )|V^| 2 ,i;), V»e% (2.7) 



r/2 

or by an explicit Euler scheme 

V2 _ .,0 



7 7/^^) + (Vuo, Vu) = (a(u )\V<S>° h \ 2 ,v), WveV h (2.8) 

By the classical finite element theory for elliptic equations and for interpolation, we have 

||£/° - u \\ L2 + - 4>°\\ L i2 /5 < Ch 2 . (2.9) 

Here we assume that the solution of the initial/boundary value problem fll.lR — fjl.3j) exists 
and satisfies 

||?Xo||^2 + \\u\\l™((0,T);H 2 ) + IKIIl°°((0,T);# 2 ) 

+ IKt|k°°((0,T);.H"l) + \\ u u\\l 2 ((0,T);H 2 ) + \\ u ttt \\l 2 ({0,T);L 2 ) < C, (2.10) 

and 



\L°°({0,T);W 2 < 12 / 5 ) + II < ^I|l°°((0,T);W 1 > 6 ) + I|V0||l°°((O,T);L°°) 

+ \\9\\L°°((o,Ty,w 2 - 12 / 5 ) + \\^9\\L°°((o,Ty,L°°) ^ C- ( 2 -ll) 

The emphasis of this paper is on the unconditionally optimal error analysis. The above 
regularity assumptions may possibly be weakened for the analysis below. We present our 
main result in the following theorem. The proof will be given in Sections 3-4. 
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Theorem 2.1 Suppose that the system il.l\) - (L2]) with the initial and boundary conditions 
\1.3\) has a unique solution (u,4>) satisfying i2.10\) - f2.11]) . Then the finite element system 
\2.4\) - (2~l\) admits a unique solution (U£, <J>^ ^ 2 ), n = 1, • • • , N, such that 

max \\U% - u(-,t n )\\ L 2 + max ||<" 1/2 - </>(•, t n _ 1/2 )\\ Ll2/s < C(r 2 + h 2 ), (2.12) 

l<n<N l<n<N 

max \\U% - u(;t n )\\ m + max \\^ 1/2 - </>(;t n _ 1/2 )\\ wl , iy5 < C(r 2 + h). (2.13) 

l<n<N l<n<N 

To prove Theorem 12.11 we introduce a time-discrete system of equations: 

V • (a(U n+1/2 )V$ n+1/2 ) =0, for n = 0, 1, ■ ■ ■ (2.14) 

D T U n+1 - AU n+1/2 = <r(c> +1/2 )|W n+1/2 | 2 , forn = 0,l,--- (2.15) 

subject to the boundary/initial conditions 

U n {x) = 0, $ n+1/2 (x) = g(x,t n+l/2 ) for x £ 30, 1 < n < N, 

U°(x) = uq(x) for x £ f2 

771/2 _ 770 _ (2.16) 

— At/ 1 / 2 = a(u )\V<j) \ 2 for xeSl, 

^ r / 2 

u 1 ' 2 = 0, for x e ao, 

where i^ is the solution to the elliptic equation 

f V • (ct(u )V(/> ) = in 0, 



AO, 



x) = g(x, 0) for x £ d£l. 



(2.17) 



With the solution of the time-discrete system (U n ,<& n 1 ' 2 ), we have the following error 
splitting: 



\U'h-u n \\ < \\e n \\ + \\U n -UJl\\ 

| $ n-l/2 _ 0n-l/2.i < II n-l/2| + || $ n-l/2 _ $?-l/2. 



e n = C/ n - u n , v n ' 1/2 = $ n_1/2 - <J) n ' 1/2 ■ 



where 



Note that the fully discrete system f|2.4[) - (|2.5|) can be viewed as the spatial discretization 
of the elliptic system (|2.14p - (|2.15p . The key issue is to prove the regularity of the solution 
to the time-discrete equations (|2.14p - (|2.15p required in the error estimates of the Galerkin 
finite element method. We present the estimates of the error functions (e n , n n ~ 1 ^ 2 ) and 
(U n - U%, $ n -V2 _ $™~ 1/2 ) in Section 3 and Section 4, respectively. 

Moreover, the error estimates given in the above theorem for $ are defined in the time 
level t n _i/2- To get the solution at the time level t n , we define 

n :4« +1/2 +*r 1/2 ). 



5 



By the above theorem, we see that 



max jn-H;tn)\\ L i2 /B <C(r 2 + h 2 ), (2.18) 

l<n<N 

max JM-H;tn)\\wW <C(r 2 + h). (2.19) 

Kn<N 



The following lemma can be proved by noting the definition (|2.3p and using a triangular 
inequality. 

Lemma 2.1 Let {v n }^ =0 be a sequence of functions on Q,. Then for any norm \\ ■ \\, 



T\\v n \\ <2tJ2 P m ~ 1/2 \\ + t\\v°\\ < 2VT 



m=l 



m=l 



3 Temporal error analysis 



In this section, we prove the existence and uniqueness of the solution of the time-discrete 
system (12.14p - (12.16p and establish error bounds for (e n , r/™ -1 / 2 ). 

Theorem 3.1 Suppose that the system li 1. 1]) - (T7!^) has a unique solution (u,(f>) satisfy- 
ing 12.10]) - [2.11]) . Then the time- discrete system \2.1J$ - \2.16^i admits a unique solution 
(U n , $ n - x /2) such that 

max \\U n \\ H 2+ max \\D T U n \\ H 2 + \\D T U 1/2 \\ H 2 < C, (3.1) 

l<n<N l<n<N 

max H^ 1 ^!! + max || V $n-i/2|| < c vi<p<oo, (3.2) 

l<n<7V w l<n<N 

and 

max \\e n \\ H i + max \\r) n ~ 1/2 \\ w i,i2/s < Ct 2 . (3.3) 

l<n<N l<n<N x ' 

Proof The existence and uniqueness of solution to the linear partial differential equations 
(|2.14p - (|2.16p is obvious. In the following, we only prove the estimates f)3. 1 j) - (|3.3j) . 
Since U° = u°, the error functions e n+1 and rj n+1 ^ 2 , < n < N — 1, satisfy 

- V • [a{U n+1 l 2 )Vrf +l l 2 ] = V • [(a(u n+1 ^ 2 ) - a{U n+1 l 2 ))V ^ n+1/2 ] + V • i?J +1/2 , (3.4) 

D T e n+1 - Ae n+1/2 = (a(U n+1/2 ) - a(u n+1/2 ))\V0 n+1/2 \ 2 (3.5) 
+ a(U n+1/2 )(V(j) n+1/2 + V$ n+1/2 ) • Vrf +1 / 2 + R% + \ 



and 



where 



-e 1 ' 2 - Ae 1 / 2 = i?V 2 , (3.6) 



U V2 - u ° 



Rl /2 = u t \ t=0 - + a(u°)|V«^| 2 - a(u^ 2 )\V^ 2 \ 2 

r/2 



R n+1 = fa 

U ~ dt 



D T u n+1 + A(n n+1 / 2 - u n+1/2 ) + (a(u n+1/2 ) - a(u n+1/2 ))\V(j) n+ 



1/2 12 
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R n + 1 / 2 = (a(u n+1 ' 2 ) - a{u n+1 l 2 ))V4> n+1/2 

are the truncation errors. With the regularity given in (|2.10p - (|2.1ip . we have the following 
estimates for the truncation errors: 

/N-i \ h 

o< max i lK 1/2 || L a+ 1^t\\R:+% 2 ) +r\\R]/ 2 \\ L6 <Cr 2 . (3.7) 

- U - \n=0 ) 

To prove (|3.1|) - (|3.3|) . first we study the error e 1 / 2 . Multiplying the equation (|3.6|) by 
|g-i/2|4^i/2 anc j integrating it over f2, we get 

\\e l l 2 \\ L ,<Cr 2 

which further shows that T\\e l / 2 \\ H 2 < Ct 2 . Since H 2 ^ C a , it follows that U 1 / 2 E C a (Tt). 
By applying Schauder's W x,v estimate (HE! [27] to ^ ne equation (13.4)) with n = 0, we derive 
that 

\\v 1/2 \\w^ < CWe 1 ' 2 ^ + C\\Rl /2 \\ L e < Ct 2 . 
By noting the fact e° = 0, from the equation (|3.5[) with n = 0, we get 

HVe 1 !!!, + rllAe 1 !!!, < r{\\e l l 2 \\ 2 L2 + ||(V^ +1 / 2 + V$ n+1 / 2 )V V 1/2 \\l* + H^HlO < Cr\ 
To conclude, we have 

lle 1 !!^ + |^ 1/2 ||vki.6 +r 1/2 ||e 1 ||^ 2 < C r 2 . (3.8) 



Secondly, we present L 2 error estimates for the solution of (13.4[) - (]3. 5H . Multiplying (I3.4p 
by rf l+l / 2 and integrating the result over Q, we obtain 

\\v n+1/2 \\m < C||e n+1/2 || L 2 + Ct 2 . (3.9) 
Again, multiplying (|3.5p by e™ +1 / 2 and integrating it over Q give 

iz? r (|| e " +1 ||| 2 )+||Ve" +1 / 2 ||| 2 

< C\\^ +1 / 2 \\ L 4e n+1 / 2 \\ L 4Vr +1/2 \\l* + IIK +1 HL 2 ||e n+1/2 || L2 
+ (o{U n+l l 2 )(V<t> n+l/2 + V$ n+1 / 2 )e n+1/2 , V77 n+1/2 ). 

Using (I2.15P and integrating by parts, 

\(a(U n+1 / 2 )(V<p n+1/2 + V$ n+1 / 2 )e n+1 / 2 , Vr? n+1 / 2 )| 

< \ (a(U n+1 / 2 )e n+1/2 Vcp n+l/2 , Vt/ 1+1/2 )| 

+ |(V • (a{U n+l l 2 )V^ n+1 ' 2 )e n+l/2 + a{U n+l/2 )V<5> n+1/2 ■ Ve n+1/2 , rj n+1 / 2 )\ 
= |(a(C/ n+1 / 2 )e n+1 / 2 V^ +1 / 2 , Vrf +1 / 2 )| 

+ \{a{U n+1 ' 2 )Vif +l l 2 ■ Ve" +1 / 2 , rf + 1 / 2 ) + (a{U n+1 l 2 )Vcp n+1/2 ■ Ve n+1 / 2 , t^ 1 / 2 )] 

< C'||e n+1 / 2 || L2 ||V// n+1 / 2 || £ 2 +C||Ve n+1 / 2 || L2 (||Vr ? " +1 / 2 || i 2||? 7 ™ +1 / 2 || L o + ||rf +1 / 2 || i2 ) 
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Applying the maximum principle to the elliptic equation (|2.14p . we obtain ||3 >n+1/ ' 2 ||i y °° < C 
and so ||r? n+1/2 ||L°° < C for 1 < ra < iV - 1. It follows that 

^.(l|e" +1 |li 2 )+i||Ve^/ 2 || i2 

< C||e-+ 1 /2||2 2 + c , || -^+l/2| | 2 2 + C7 || r? ^+l/2||^ i + ^H^+ln^^ 

<C(||e"+ 1 ||l 2 + ||e n ||| 2 + || e ' 1 - 1 ||i 2 )+Cr 4 

where we have noted (|3.9p and used the inequality || ■ ||^6 < C\ \ ■ \\ H i . By applying Gronwall's 
inequality to the above inequality, with (|3.9j) we see that there exists a positive constant 
t\ > such that when r < n, we have 



JV-l 

max ||e" +1 || 2 2 + max \\v n+1/2 \\m + T \\e n+1/2 \\ 2 m r < Cr 4 . (3.10) 

Kn<7V" Kn<7V— 1 MjH - v ' 



n=0 

Finally, we study the regularity in (|3.1|) - ([3.2p and the estimate for ||r7 n+1 || l/[/ i,i 2 /5. Note 
that the above estimate implies that 

n 

\\D T e n+1 \\ 2 L2 + r\\D T e m+1 ' 2 \\ 2 m < CV 2 . ( 3 - n ) 

m=0 

||D T [/ n+1 || L2 < \\D T u n+l \\ L 2 + \\D T e n+1 \\ L 2 < C 

for < n < N — 1. Regarding ()3.5[) as an elliptic equation and applying the H 2 estimate 
[TT] . with (|3.1ip we obtain 

||e n+1 / 2 |b 2 < C (\\D r e n+l \\ L 2 + ||e- +1 / 2 || L2 + llVr?^ 2 ^ + \\V V n+l/2 \\% + \\K +1 \\l*) 
<Ct + CIIV^+^II^IIv^+^h^ 

<C 1 T(l + ||V ?? n+1 / 2 ||^ / 6 2 ), (3.12) 

forn = l,--- ,N-1. 

Now we prove a primary estimate 

||V? ? n+1/2 || L 6 < 1 (3.13) 

for < n < N — 1, by mathematical induction. It is easy to see from (|3.8p that (I3.13P holds 
for n = Oifr<l/Co. We assume that (|3.13j) holds for < n < k. Then from f|3. 12j) and 
dHH]) we get 

\\e n+1/2 \\ H 2 < (2d + C r)r, for < n < k 

and by Lemma 2.1, 

n 

\\e n+1 \\ H2 < 2 l|e m+1/2 ||//2 + ||e°||if2 < (2Ci + C t)T, for 1 < n < k. 

m=0 

Hence, 

||^ n+1 ||i/2 < \\u n+1 \\ H 3 + (2Ci +C r)T < C 3 , for 1 < n < k. 
Since # 2 ^ C a in M d (d = 2, 3), we have 

||{7 fc +3/ 2 || CQ < C||[/ fc+1 || H2 + C||?7 fc || H 2 < C 4 . 



S 



With the Holder regularity of a(U k+3 / 2 ), by applying the W l ' p estimate [27] to (|3.4p for 
n = k + 1 , we obtain 



||V7? fc+3 / 2 || L6 < C 5 \\(a(u k + 3 / 2 ) -a{U k+ V 2 ))V4> k+Z/2 \\ L s + C S \\R 



k+l/2. 



<C 6 ||e* +3 / 2 || L6 + C 6 r 2 

< C 7 ||e fc+3/2 || jF/ i + C 6 x 2 (3.14) 
<C 8 (\\e k+1 \\ H i + \\e k \\ H i)+C 6 T 2 

k+l 

<C 9 J2\\e m+1/2 \\m+C 6 r 2 <C l0 r, 

m=Q 

where we have used Lemma 2.1 and (|3.10p . 

By choosing r < 1/ max{Co, Cio}, we get ||V77 fc+3 / 2 || i 6 < 1 and we complete the induc- 
tion. Thus, we have proved that (|3.13p holds for < n < N — 1, which together with (|3.12p 
implies that maxi< n <Ar„i ||e n+1//2 ||^2 < Cr. By using Lemma 2.1 again, we obtain 

max \\U n \\ H i <C. (3.15) 

l<n<Af 

Since H 2 <^-» C a in M d (d = 2,3), with the Holder continuity of [/ n + 1 / 2 ) we apply the 
W 1 '? estimate [27] to ([2TT4"]) and derive that 

max ||V$ n+1/2 || L P < C p , Vl<p<oo (3.16) 

where we have noted H^H^/i.p < C. With the estimates (|3.15p - (|3.16p . we can perform the 
W 2 ' p estimate (with p = 12/5) [HE] for the elliptic equation (|2.14p to obtain 

max ||$ n+1 / 2 \\ W 2,i2/s < C (3.17) 

l<n<N— 1 

where we have also noted ||<7|| ^2,12/5 < C. 

With the above estimates, multiplying (|3.5p by — Ae n+1 / 2 and using the inequality 
||gn+i/2||^ 2 < (7|| Ae n+1//2 || L 2 (because of the boundary condition e n+l l 2 = on dQ), we 
obtain 

Dr\\e n+1 \\ 2 m + \\e n+1/2 \\ 2 H^ 

< C(\\Vr +1/2 \\l~ + l|V* n+1/2 ||i3)||Vr, n+1 / 2 ||2 6 +C \\e n + 1 \\ 2 L2 + C \\K +1 \\ 2 L2 
<C\\^ +1 I 2 \\ 2 h1 +Ct\ 



where we have used (|3.14|) . By Gronwall's inequality, we see that 

N-l 

lie 11 jji -t- 

0<n<N-l 



max ||e n+1 || 2 ffl + y r||e n+1 / 2 || 2 2 < Cr 4 (3.18) 

n=0 

and by using (I3.14D again, we have 

|| Vr? n+1 / 2 || il2/5 <Ct 2 . (3.19) 
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Moreover, by (|3.18j) and Lemma 12. 11 we have further 

\\e^ H ,<Cr^(^r\\e k+ yY H X 2 <Cr, 

which implies that ||-D r e n ||^2 < Cr _1 ||e n ||^2 < C. With the regularity assumption (|2,10p . 
we see that 

\\D T U n \\ H 2 < \\D T e n \\ H 2 + \\D T u n \\ H 2 < C. (3.20) 

So far we have proved that there exists a positive constant To such that for r < tq 
(|3.ip - (|3.3p hold. Also we have proved that for any r > 0, (|3.ip - (|3.3p hold at the initial 
steps. 

For r > to, if we assume that 

WV^-^Wlp + \\U n \\ H 2 < C n (3.21) 

for 1 < n < k (for any fixed 1 < p < oo), where C n is a constant dependent upon n. Then, 
by writing the equation (|2.14p as 

(l - ^AjU k+1 = (l + -&)U k + <?{U k+l,2 )\V§ k+l,2 \ 2 T, (3.22) 

and applying the classical W 1,p and W 2 ' p estimates [II 111} [27] of elliptic equations to (|2,15p 
and (1332]) . we get 

||V$ fc+1 / 2 || LP + ||C/ fe+1 || H2 <C fe+1 , 

where Ct+i depends upon C n , 1 < n < k. By mathematical induction, (|3.2ip holds for 
1 < n < N. Since N < [T/tq] + 1, by setting C* = maxi< n <7v{C n } we obtain 

|| V$ n+i/2|| Lp + ||[/«+i||^ 2 < c*, n = 1,2,- ■■ ,N. 
(|3.1j )-( j3~3j) follow immediately. The proof of Theorem 13.11 is complete. I 



4 Spatial error analysis 

In this section, we present error estimates of the Galerkin finite element method for the 
time-discrete system (l2TT4l - (l2~T5l) . Let = g° + U h (<p° - g°) and p™ -1 / 2 ^-!^ = 

g n-l/2 + n h (*n-i/2 _ g" 1 - 1 / 2 ) for n = 1, 2, • • • , N, and define R h : flj(fi) V h to be a 
Riesz projection operator defined by 

(V(w - R h v), Vw) = 0, for all v E flj(fi) and w G V h . 

We summarize some basic inequalities below. The proof follows the classical finite element 
theory for elliptic equations, see [HI [33] for references. 

\\w\\ W m,p <Ch wp - d/q *>\\w\\ W m, q , w G Vh, l<q<p<oo,m = 0,l, (4.1) 
IIHIw 1 -? — C 1 1 1 w l I ; 1 < P < oo, (4-2) 
\\v -U h v\\ LP + h\\v -U h v\\ w i, P < Ch 2 \v\ W 2, P , p>d/2 (4.3) 
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< Ch 2 ^ 1 - 1 ' 2 - g n -V 2 \ W 2, P , p > d/2, (4.4) 
\VU h v\\ LP <C\\v\\ w i, P , for all ^^(fl) with p>d, (4.5) 



and 



\\Rhv\\ w i, P < C\\v\\ w i, P , for all v E W l * ', 1< p < oo, (4.6) 

\\v - R h v\\LP + h\\v - R h v\\ w i, P < Ch 2 \\v\\ W 2, P , for all v£\V 2 ' p , 1 < p < oo, (4.7) 

v ||ty2, g , dp/(d + 2p) < q < p. (4.8) 

Let r/° = $° - P°c/> and 

e h l = U%-R h U n , = $ «-i/2 _ p n-i/2 $n _i/ 2; for 1 < n < iV _ 

We present error estimates of the spatial discretization in the following theorem. 

Theorem 4.1 Suppose that the system M. 1\) - (T73\) has a unique solution (u,(j)) satisfying 

\2.1U\) - F2,.llS) . Then the fully- discrete finite element system ^2.4\ )- [2~7\ ) admits a unique 
solution (U'h, ), n = 1, 2, • • • , N , such that 



max hhWi? + max ||% 1/2 H^ia/s < C/i 2 , (4.9) 

l<n<iv l<ri<iV 

n-1/2. 



max ||Ve^|| i2 + max ||V^ ^\\ L12/ s < C/i. (4.10) 

\<n<N l<n<N 

Proof At each time step of the scheme, one only needs to solve two uncoupled linear discrete 
elliptic systems. It is easy to see that coefficient matrices in both systems are symmetric 
and positive definite. The existence and uniqueness of the Galerkin finite element solution 
follows immediately. Since the inequality (|4.10p follows from (|4.9p via the inverse inequality 
AO) , it suffices to prove (|4Uj) . 

Let $° = The solution of the time-discrete equations (|2.14|) - (|2.15|) satisfies 

{D T U n+1 , v) + (V[T +1/2 , Vv) = {a(U n+1 / 2 )\V^ n+1 / 2 \ 2 , v), (4.11) 
(a(U n+1 / 2 )V<S> n+1 / 2 , V(p)=0, n = 0, 1, • • • , N (4.12) 

for any v, <p G Vh, and 

( ^p , v) + (W 1 ^ Vv) = (aK)|V/| 2 , v), vGV h (4.13) 

From the above equations and the corresponding finite element system (|2.4j) - (|2.7p . we find 
that the error functions e^ +1 , r/^ 1 ^ 2 £ Vh, < n < N, satisfy 

(D T el+\ v) + (Ve^ 2 , V.) 

= (D T (U n+1 - R h U n+1 ), v) + {(o-(U^ +1/2 ) - a(U n+1 ^ 2 ))\W +1 / 2 \ 2 , v) 
+ 2((a(UZ +1/2 ) - o-{U n+l / 2 ))V<$> n+1 l 2 ■ V(^ +1/2 - Q"* 1 / 2 ), v) 
+ Hl); +1/2 )|V« +1/2 - $«+V2 ) |2 > w) 
+ 2(a(f/ n+1 / 2 )V$ n+1 / 2 • V($™ +1/2 - $™ +1 /2) 5 „) 
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J2ir l/ \v), (4-14) 



i=l 



(a(C/™ +1 / 2 )V< +1/2 , Vy>) = -((a(U^ +1/2 ) - a(U n+1 ' 2 ))V^ 1/2 , Vp) 

+ (<7(£>+ 1 / 2 )V(<I> n + 1 / 2 - P^+Va^n+i^^ V(p j ; (4 15) 



and 



-1/2 



r/2 ' 

(a( Uo )(|V^°| 2 - \V<S>° h \ 2 ),v) + (l^C/ 1 / 2 - R h D T/2 U l '\v), (4.16) 



for all u, (f € V^, where Z?^?/ 1 / 2 = 2(C/ 1 / 2 - u°)/r. 

First, we estimate the error functions at the initial step. Since V • (cr(u )V</> ) = in Q,, 
by setting v = e ft in (14.16P we have 



D T/2 fr' 2 - R h D T/2 U l /\e 1 ' 2 ) < C\\D T/2 U^\\ H ,h 2 \\^ 2 \\ L , < e||e£ /2 || H i + Ce~ l h A 



and 

|(aM(|V0°| 2 -|V$°J 2 ),gi /2 )| 

< |(o-(«o)|V(*S - /)| 2 , ef )| + 2|(^ )V<1> • V($° - /), ef )| 
= |(a(n )|V($° - /)| 2 , ef )| + 2| (a(u )V$°($° - 0°), Vef )| 

< C||eJ/ 2 || L6 ||V(/ - $°)|| 2 12/5 + C||^ - /|| L2 ||Vef || L2 

< e||ei /2 || H i +Ce~ 1 /i 4 , 

where we have used integration by parts and (|2.9|) . With the above estimates, (|4.16|) reduces 
to 

\\t' 2 \\ 2 L2 <\\el\\l 2 +Crh'<Ch'. (4.17) 

Secondly, we present estimates for ||e^ +1 ||£ 2 and ||^ +1//2 ||ii2/5 for < n < N — 1. For 
this purpose, we take v = e^ +1 ^ 2 in (|4.14p and we have 

/ n + l/2 ( _n + l/2 ) < lle^ll^n^^+l.^^n+lll^ 

< C||V^ +1/2 || L2 ||D r C/ m+1 - i^t/" 4 " 1 ^ 
^ellVe^lH.+Ce- 1 !!^^ 1 !! 2 ,^ 4 , 

/ n + l/2 ( _n + l/ 2) < C( ||^+l/2|| i2 + ||£n + l/2 - ^[/^ j,^) || V $"+V2||2 g ||^ + l/2 

<C(||^ +1/2 || i2 + / l 2 )||Ver 1/2 || L2 



< e ||Ver 1 / 2 ||i 2 +C e - 1 (|lC 1/2 Hi 2 +^ 4 ) ) 

/ n. + l/2 ( _n + l/2 ) < c|| _n+l/2 ^ ( ^ + n+1/ 2 _ (4 1§) 

• (l|V ?? r 1/2 |l^ + ||V(d»" +1 /2 _ p; +1 /2^ + l/2 )||L6)||v<J) n + l/2 ||L6 
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< ^||Ver 1/2 ||L 2 (llC 1/2 |lL 2 + ^)(^ /4 ||Vi 1+1/2 || il2/5 + W +1/2 - 9 n+1/2 \\ H >) 

< e ||Ve- +1/2 ||| 2 + Ce-Hll^ 172 !!!, + ^X^-^ll^r 172 !!^^ + c), 
i2 +1/2 (4 +1/2 ) < c||^ +1 / 2 iu e ||v($" +1 /2 _ < +1 / 2)|| 2 i2/5 

<C||Ver 1 / 2 || L2 (||Vr ? r 1 / 2 ||2 12/5+ / l 2 ) 

< e iiver 1/2 ii! 2 +^" i (iivr ? r i/2 ni 12/s +/ l 4 ) 

/ n + i/ 2( _n + i/ 2) < ]K +1/2 \\ L 4K +1/2 ~ * n+1/ X^4^ n+1/2 \\ L ^ 



< e||V^ +1/2 ||| 2 + Ce- 1 ^ 1 ' 2 - ^ +1 /2|| 



2 

£ 12/5 • 



where we have used (|4.5|) in the estimate of /3 +1//2 (e^ +1 ^ 2 ). 

By applying the W 1 '^ estimate [26] to (|4.15p and using (|4.4p . we obtain 

l|V< +1/2 ||£i 2/5 (4.19) 

< C (||(<t(C^ +1/2 ) - a(U n+1 / 2 ))V$ n h +1/2 \\ Ll2/5 + ||a(^ +1 /2 )V(<i) n + i/2 _ p n+l/2^ n+1/2)hi2/B 

< C||C/™ +1/2 - f> n+1/2 ||£4(||V^ +1/2 || L6 + ||VP, n+1/2 ^ +1 /2|| L6 ) + c/* 

< C(\\r h +1/2 \\ Li + ||[/" +1 / 2 _ i^+V2|| L4)( || V ^+l/2 ||L6 + || V pn+l/2 $ n + l/2 ||L6) + ^ 

< Ch- d '\\\^ 2 \\ L 2 + / i 2 )(/ i - d/4 ||V ?? " +1/2 || L12/5 + C) + C7i 

< c^ii^ii^nv^iLiva + h-^\\r h +1 ' 2 \\ L , + h^\\vC ll X^ + h) , 

where we have used the inverse inequality (|4.ip and (|4.8|) with q = 2 and p = 4. To estimate 
||7/ n+1 / 2 || L i2/5, we rewrite (|4.15p as 

(a(f> n+1 / 2 )V(^ +1 / 2 - ^ +1/2 ), V<^) + (Hl> +1 / 2 ) - a(UZ +1/2 ))V<S> n h +1/2 , V<p) = 
and apply the Nitsche technique. Define ip as the solution of the elliptic equation 

-V • {a{U n+l / 2 )V^) = |$ n+1 / 2 - ^ +1 /2|2/5 ($ n+l/2 _ $ n+l/2 ) 

with the boundary condition tp = on 30. The solution ^ to the above elliptic equation 
satisfies that 

U\\ W 2, l2/ r < C\\^ 2 - ^ 2 \\f 2/5 . 

Since 



, $ n+l/2 _ 


^n+l/2,,12/5 
II £18/5 










= {o0 n - 


-V2)V($"+V2 


_ $ n+l/ 2) 




= {°(u n - 


-V2)V($«+V2 


- < +l/2 ) 


, V(^ 


= {°(u n - 




_ $ n+l/ 2) 


, V(^- 



fn+1/2% _ ^^^+1/2^^7^+1/2 



< c||v(^ +1 / 2 - ^ +1/2 )||£ 12/5 ||v(^ - n^)|| L12/7 

+ C ||^+l/2_^+l/2 ||L2||w n+l/2 ||L4||vn ^ ||L4 

<C/ t ||V(^ +1 / 2 -^ +1 / 2 )|| il2/5 ||^|| w2 , 12/7 

+ ||^/2 _ ^+1/2 ||l2 / ||v?? n+l/2 ||L4 + l]V pn + y^l/2 h \ 
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|v(^ +1 / 2 -< +1/2 )|| L12/5 



< C\\tp\\ W 2,12/7 

_i_ ||r>n+l/2 fT n+l l 2 \\ (h ~ 7 II V7^ n+1 / 2 II _l II rT>™+l/2 11 , 11 n+l/2|| \ 

+ 11^ ~ ^/l III, 2 MlV^ || L 12/5 + ||* ' || W 2,12/5 + ||5 H^/2,12/5 1 

we derive that 

\\$ n+1/2 -$ n h +1/2 \\ L iy S (4.20) 
< Ch\\V(<5> n+l / 2 - < +1/2 )|| L i 2 /5 + C||l> +1 / 2 - UZ +1/2 \\ L 2 (/i- d/6 ||V^ +1/2 || Ll2/5 + C) . 

By flUIlD and (|4TTg ]) - (|430j) . (jlTlD reduces to 

^(llerili^ + HVer 1 / 2 !!^ (4.21) 
< Ce^\\D T U^f Hl ^ + C (\\e n h \\l> + hVWh) 

+ c (^/ 3 ||v< +1 / 2 ||^ 12/5 ||^v 2 ||2 2 + ||v^ +1/2 lll 12 /. + ^ 4 ) 

forn = l,--- ,N-1. 

Now we prove a primary estimate 

\\el\\ L 2 <h 7/4 , 0<n<N (4.22) 

by using mathematical induction. 

By (|2.9p and (|4.7p . this inequality holds for n = if h < h\ for some given positive 
constant h\. If we assume that (|4.22p holds for < n < k, then from (|4.17p we know that 
||e^ +1 ^ 2 || i 2 < 2/i 7 / 4 + Ch 2 for < n < k and from the inequalities (|4.19p we see that there 
exists a positive constant /12 such that when h < h,2, 

l|V^ +1/2 || L i2/ 5 < Ch, for 0<n<k. (4.23) 
With the above inequalities, (|4.2ip reduces to 

Dr(K +1 \\h) + HVer 1/2 ||i 2 < C{\\el\\l> + \\ef ||| 2 ) + Ch* 
for < n < k, which implies that there exists a positive constant C\\ satisfying 

\K +1 \\ L 2 < C u h 2 (4.24) 

for < n < k. Therefore, ||e^ +1 ||^2 < /1 7 / 4 if ft, < mm.(h\, /12, 1/Cf-^), which completes the 
induction. Thus, (02^- ([EM]) also hold for all < n < N - 1, and from (|4"~20"|) we derive 
that ||$ n+1 /2 - $^ +1/2 || L i2/5 < Ch 2 and so 

||^ +1/2 || L1 2/s < ||^ +1/2 - ^ +1/2 Hl-/ 5 + ||P; +1/2 ^ +1 / 2 - ^ +1 / 2 || Ll2/5 < Cft 2 (4.25) 
for < n < N - 1. 

So far we have proved that the estimate (|4.9p holds if h < /io> f° r some positive constant 
ho. It remains to show that 

max ||e^|| L2 + max ||^ +1/2 || Ll2/5 < C (4.26) 

0<n<N 0<n<N—l 
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for h > ho. In fact, substituting tp = c &^ +1 ^ 2 in (j2.4j) we get 

||V<I>™ +1/2 || L2 <C, n = l,2,---,N. 
By the inverse inequalities, we have 

l|V$£ +1/2 |k°° < Chf, n = l,2,---,iV 

and so 

l|V< +1/2 || L i 2 /s <C fto , n = l,2,--- ,iV. 
where C^ is a positive constant dependent upon Iiq. From (|2,4|) with u = f/ n+1 / 2 we g e t 

^(iic/nno + \m n h +1/2 \\h < cuv^n^irr 172 !!^ 



^ciiv^m. + Invur 1 / 2 !!!,, 

which implies that 

||t/r 1 |l! 2 <||^|l! 2 +C7, r<...<C, o r. 
This completes the proof of Theorem 14.11 I 

Theorem 12.11 follows from Theorem 13.11 Theorem 14.11 together with (|4.4p and (|4.7p . I 

5 Numerical results 

In this section, we present two numerical examples to illustrate our theoretical results. The 
computations are performed with the software FEnics. 

Example 4-1- We rewrite the system (|l,ip -( fL2j) by 

^-Au = a(u)\V<j ) \ 2 + f 1 , (5.1) 
-V • (<7(u)V0) = h, (5.2) 

where O = (0, 1) x (0, 1) and 

1 + u A 

The functions /i, /2, and the Dirichlet boundary conditions are chosen corresponding to 
the exact solution 

u(x, y, t) = exp(x + y — t), 4>(x, y, t) = 1 + sin(x + y + t). 

A uniform triangular partition with M + 1 nodes in each direction is used in our compu- 
tation. We solve the system by the linearized Crank-Nicolson Galerkin method with linear 
elements and quadratic elements, respectively. To confirm our error estimates in the L 2 
norm, we choose r = h for the linear FEM and r = /i 3 / 2 for the quadratic FEM. We present 
the numerical results in Tables [T][2j We can see clearly from Tables [U that the L 2 errors of 
the linear FEM are proportional to h 2 and from Table [2] that the L 2 errors of the quadratic 
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FEM are proportional to h 3 . To see the errors in the H 1 norm, we take r = /i 1 / 2 for the 
linear FEM and r = /i 3 / 2 for the quadratic FEM, and we present numerical results in Table 
OH All these results are in good agreement with our theoretical analysis. 

To show the unconditional stability, we test the linearized Crank-Nicolson Galerkin 
method with linear elements, h = 1/80 and the large time steps r = h, 5h, lOh, 20h. We 
present numerical results in Table [5j The results show that the scheme is stable for large 
time steps, although the numerical results with r = 20h seem not very accurate. 

Table 1: L 2 errors of linear FEM with h = r = l/M (Example 4.1). 







\\U%-u(;t n ) 








t 


M = 20 


M = 40 




M = 80 


order 


1.0 


7.8063e-05 


1.9587e-05 




4.9042e-06 


2.00 


2.0 


9.9117c-05 


2.4975e-05 




6.2605e-06 


1.99 


3.0 


8.7998e-05 


2.2134e-05 




5.5422e-06 


1.99 


4.0 


5.6591e-05 


1.4242e-05 




3.5666e-06 


1.99 














t 


M = 20 


M = 40 




M = 80 


order 


1.0 


7.2691e-05 


1.7791e-05 




4.3746e-06 


2.03 


2.0 


9.7524e-05 


2.3836e-05 




5.8610e-06 


2.03 


3.0 


1.3954e-04 


3.4376e-05 




8.4930e-06 


2.02 


4.0 


1.4342e-04 


3.5393e-05 




8.7511e-06 


2.02 



Table 2: L 2 errors of quadratic FEM with h = l/M and r = /i 3 / 2 (Example 4.1). 







\K-u(;t n ) 








t 


M = 10 


M = 20 




M = 40 


order 


1.0 


8.8214e-06 


1.2113e-06 




1.5496e-07 


2.92 


2.0 


2.0196e-05 


2.5321e-06 




3.1437e-07 


3.00 


3.0 


1.7212e-05 


2.1149e-06 




2.5866e-07 


3.03 


4.0 


4.7866e-06 


5.5294e-07 




6.3392e-08 


3.12 






\\n-<t>(;t n )\ 








t 


M = 10 


M = 20 




M = 40 


order 


1.0 


3.4542e-05 


4.0535e-06 




5.0059e-07 


3.05 


2.0 


3.9065e-05 


4.7335e-06 




4.3953e-07 


3.24 


3.0 


4.1164e-05 


3.2571e-06 




6.2613e-07 


3.02 


4.0 


3.0560e-06 


2.5956e-05 




3.7303e-07 


3.06 
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Table 3: H 1 errors of linear FEM with h = l/M and r = h 1 / 2 (Example 4.1). 





\\UZ-u(;t n )\\ H i 


t 


M = 40 


M = 80 


M = 160 


order 


1.0 


5.6024e-03 


2.3706e-03 


9.9037e-04 


1.25 


2.0 


3.6159e-03 


1.8179e-03 


7.7195e-04 


1.11 


o.u 


2.8675e-03 


1.3254e-03 


6.4766e-04 


1.07 


4.0 


1.7923e-03 


8.2924e-04 


3.5460e-04 


1.17 




\\n-^-^n)\\m 


t 


M = 40 


M = 80 


M = 160 


order 


1.0 


4.1118e-03 


2.2510e-03 


1.1312e-03 


0.93 


2.0 


4.8396e-03 


2.0593e-03 


1.3143e-03 


0.94 


3.0 


5.0531e-03 


2.7701e-03 


1.0180e-03 


1.16 


4.0 


5.1689e-03 


1.9872e-03 


1.4263e-03 


0.93 



Table 4: H 1 errors of quadratic FEM with h = r = 1/M (Example 4.1). 





\\UZ- U (;t n )\\ H l 


t 


M = 10 


M = 20 


M = 40 


order 


1.0 


1.6700e-03 


3.1162e-04 


5.9398e-05 


2.41 


2.0 


1.3279e-03 


2.8247e-04 


6.4458e-05 


2.18 


3.0 


1.0156e-03 


2.2534e-04 


5.2879e-05 


2.13 


4.0 


5.0056e-04 


9.5814e-05 


1.9466e-05 


2.34 




\\n-^^n)\\m 


t 


M = 10 


M = 20 


M = 40 


order 


1.0 


1.8020e-03 


4.4199e-04 


1.0614e-04 


2.04 


2.0 


1.6539e-03 


4.0090e-04 


9.4222e-05 


2.07 


3.0 


1.6141e-03 


3.8446e-04 


8.8813e-05 


2.09 


4.0 


1.5725e-03 


3.7110e-04 


8.5031e-05 


2.10 



Table 5: L 2 Errors of linear FEM with h = 1/M and r = kh (Example 4.1). 





\\UZ- U (;t n )\\ L l 


t 


k = l 


k = 5 


k = 10 


k = 20 


1.0 


4.9042e-06 


3.6772e-05 


2.2480e-04 


1.6693e-03 


2.0 


6.2605e-06 


8.0524e-05 


3.1973e-04 


1.5965e-03 


3.0 


5.5422e-06 


6.6692e-05 


2.6145e-04 


1.1952e-03 


4.0 


3.5666e-06 


1.7407e-05 


6.6093e-05 


3.9786e-04 




\\n -0Mn)||L2 


t 


k = 1 


k = 5 


k = 10 


k = 20 


1.0 


4.3746e-06 


1.4126e-04 


5.5525e-04 


1.3714e-03 


2.0 


5.8610e-06 


1.2561e-04 


4.9591e-04 


1.0407e-03 


3.0 


8.4930e-06 


1.1959e-04 


4.6941e-04 


8.5232e-04 


4.0 


8.7511e-06 


1.1356e-04 


4.4541e-04 


7.2633e-04 
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Example 4. 2. In the second example, we consider the system (|5. 1|) - (|5.2|) in the three- 
dimensional space with the exact solution 

u(x, y, z, t) = exp(2x + y — z)(2t + sin(t)), (5-3) 
4>(x, y, z, t) = sin(x — 2y) cos(z) exp(i) . (5-4) 

We use a uniform tetrahedral partition with M + 1 nodes in each direction (see Figure 
1). The total number of tetrahedra is 6M 3 and the total number of vertices is (M + 1) 3 . We 
solve the system by the proposed Crank-Nicolson Galerkin method with linear elements. 
Table E] contains the 1? errors of the numerical solution with r = h and h = 1/10, 1/20, 1/40. 
Similarly, we can see that the I? errors for both u and </> are proportional to h 2 . 

Previous analysis for the three-dimensional problem often requires a stronger time-step 
condition than that for the two-dimensional problem. Finally, we test the linear Galerkin 
method with h = 1/80 and large time steps r = h,5h, lOh. The results are presented in 
Table [71 Numerical results show that the scheme is unconditionally stable. 




Figure 1: The three-dimensional mesh (Example 4.2). 
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Table 6: I? Errors of linear FEM with h = r = l/M (Example 4.2). 





\\U%-u{;t n )\\ L . 


t 


M = 10 


M = 20 


M = 40 


order 


1.0 


1.1089e-03 


2.8319e-04 


7.1320e-05 


1.9793 


2.0 


8.6316e-04 


2.2523e-04 


5.6987e-05 


1.9605 


3.0 


4.0520e-04 


1.0626e-04 


2.6895e-05 


1.9566 


4.0 


3.6125e-04 


9.4243e-05 


2.3822e-05 


1.9613 






t 


M = 10 


M = 20 


M = 40 


order 


1.0 


4.0129e-04 


1.1562e-04 


2.9779e-05 


1.8761 


2.0 


7.8577e-04 


2.1723e-04 


5.5611e-05 


1.9103 


3.0 


7.7231e-04 


2.1482e-04 


5.5087e-05 


1.9047 


4.0 


5.0533e-04 


1.4801e-04 


3.8404e-05 


1.8589 



Table 7: 1? Errors of linear FEM with h = 1/M and r = kh (Example 4.2) 





\\U%-u(;t n )\\ L2 


t 


k = l 


k = 5 


k = 10 


1.0 


7.1320e-05 


2.5286e-04 


2.3962e-03 


2.0 


5.6987e-05 


2.4133e-04 


1.6463e-03 


3.0 


2.6895e-05 


1.0359e-04 


8.5356e-04 


4.0 


2.3822e-05 


2.6147e-04 


1.1001e-03 




\\n-<i>(->tn)\\v 


t 


k = 1 


k = 5 


k = 10 


1.0 


2.9779e-05 


2.7311e-04 


4.6730e-04 


2.0 


5.5611e-05 


1.2094e-04 


1.3614e-03 


3.0 


5.5087e-05 


1.0514e-04 


1.7040e-03 


4.0 


3.8404e-05 


1.0333e-04 


1.6532e-03 



6 Conclusions 

We have presented an uncoupled and linearized Crank-Nicolson Galerkin finite element 
method for the nonlinear time-dependent thermistor equations in the ci-dimensional space 
[d = 2, 3) and provided unconditionally optimal error estimates in both I? and H 1 norms, 
while existing analysis requires certain time-step restrictions. Our numerical results confirm 
our analysis and show that the proposed scheme is efficient. Our approach presented in this 
paper can be extended to many other nonlinear parabolic systems, high-order finite element 
approximations and other time discretization schemes, while the analysis only focuses on 
the electric heating model with a linear finite element method to illustrate our idea. 
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